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Abstract. Many cellular and subcellular biological processes can be described in 
terms of diffusing and chemically reacting species (e.g. enzymes). Such reaction- 
diffusion processes can be mathematically modelled using either deterministic partial- 
differential equations or stochastic simulation algorithms. The latter provide a more 
detailed and precise picture, and several stochastic simulation algorithms have been 
proposed in recent years. Such models typically give the same description of the 
reaction-diffusion processes far from the boundary of the simulated domain, but 
the behaviour close to a reactive boundary (e.g. a membrane with receptors) is 
unfortunately model-dependent. In this paper, we study four different approaches 
to stochastic modelling of reaction-diffusion problems and show the correct choice of 
the boundary condition for each model. The reactive boundary is treated as partially 
reflective, which means that some molecules hitting the boundary are adsorbed (e.g. 
bound to the receptor) and some molecules are reflected. The probability that the 
molecule is adsorbed rather than reflected depends on the reactivity of the boundary 
(e.g. on the rate constant of the adsorbing chemical reaction and on the number of 
available receptors), and on the stochastic model used. This dependence is derived for 
each model. 

Keywords: stochastic simulation, boundary conditions, reaction-diffusion problems. 
1. Introduction 

Let us consider a system of k chemicals diffusing and reacting in a domain DcK 3 . Let 
nj(x, t), j = 1, . . . , k, be the density of molecules of the j-th chemical species at the 
point x£fl. Assuming that there are a lot of molecules present in the system, the time 
evolution of density rij (x, t) can be computed by solving the system of reaction-diffusion 
partial-differential equations 
dfi ■ 

—l = D j S7 2 n j + R j (n 1 ,n 2 ,...,n k ), j = l,...,k, (1) 

where Dj is the diffusion constant of the j-th chemical species, V 2 is the Laplace 
operator and reaction term Rj(rii, n 2 , . . . , rife) takes into account the chemical reactions 
which modify the concentration of the j-th chemical species. To describe uniquely the 
time evolution of the system, we have to introduce suitable boundary conditions for the 



Boundary conditions for stochastic simulations 



2 



system of equations (fTT). The simplest boundary conditions can be formulated in terms of 
vanishing density n^x, t) on the boundary of Q or vanishing flux through the boundary 
of Q. Coupling system of equations (CQ) with such a boundary condition, we can compute 
densities n.,-(x, t) at any time t from the initial densities n 3 (x, 0), j = 1, . . . , k. 

Reaction-diffusion processes in biology often involve low molecular abundancies of 
some chemical species. In such a case, the continuum deterministic description ([1]) is 
no longer valid and suitable stochastic models must be used instead. Various stochastic 
simulation algorithms have been proposed in the literature [U [TUJ, [Hi [21] . In general, 
the stochastic treatment of diffusion and first-order reactions (such as degradation or 
conversion) is well understood. There is less understanding (and stochastic models 
differ) when second-order chemical reactions are taken into account, e.g. reactions in 
which two molecules collide for the reaction to take place. Another important problem 
is the implementation of the correct boundary conditions for the stochastic simulation 
algorithms. On the one hand, the simple boundary conditions mentioned above are 
easy to reformulate in the stochastic case - the vanishing density on the boundary of Q 
simply means that a diffusing molecule is removed from the system whenever it hits the 
boundary; and the vanishing flux through the boundary means that a diffusing molecule 
is reflected whenever it hits the boundary. On the other hand, more realistic boundary 
conditions have to be handled with care. They can be formulated in terms of the 
partially adsorbing boundary, which means that some molecules hitting the boundary 
are reflected and some are adsorbed. The partially adsorbing boundary corresponds to 
the so-called Robin boundary condition of the macroscopic partial-differential equation 
([1]). However, this correspondence is model-dependent. 

We will see later, in Section [U that the derivation of the correct boundary 
condition depends on the stochastic model of the diffusion but not on the stochastic 
model of the chemical reactions in the solution. Consequently, we start this paper 
by studying stochastic models of diffusion only. In Section [2} we introduce four 
different stochastic approaches to model molecular diffusion and we state the appropriate 
boundary conditions. In Section [3], we present illustrative simulations of all four models, 
validating the boundary conditions presented. Moreover, we also clearly illustrate that 
the boundary conditions are indeed mo del- dependent. In Section HI we present the 
mathematical derivation of the boundary conditions for each model, i.e. we provide 
a theoretical justification of results from Section [2j Moreover, we show that all four 
models are suitable for modelling diffusion far from the reactive boundary. Section H] 
is intended for a more theoretical audience and can be skipped by a reader who is not 
interested in the mathematical justification of the boundary conditions and stochastic 
models. In Section [5j we show that reaction-diffusion models can be treated using the 
same boundary conditions which were previously derived for the corresponding models 
of the diffusion only. We conclude with summary and outlook in Section [6j 
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2. Boundary conditions for stochastic models of diffusion 

The boundary condition of any stochastic simulation algorithm can be formulated as 
follows: whenever a molecule hits the boundary, it is adsorbed with some probability, 
and reflected otherwise. The special cases of this boundary condition are: (a) the 
molecule is always reflected (such a boundary is called the reflecting boundary in what 
follows); and (b) the molecule is always adsorbed (in this case the boundary is called 
fully adsorbing). The reflecting boundary condition is often used when no adsorption of 
the diffusing molecules on the boundary takes place. On the other hand, if the molecule 
can chemically or physically attach to the boundary, then one has to assume that the 
boundary is (at least) partially adsorbing. 

The important question is, what is the probability that the particle is adsorbed 
rather than reflected, and how does this relate to the reactive properties of the boundary 
for a given stochastic model? To answer this question, let us follow the x-coordinate of 
the diffusing molecule (the other coordinates can be treated similarly), so that we study 
the diffusion of molecules in the one- dimensional interval [0, L\ where L is the length of 
the computational domain. Assuming that we have a lot of molecules in the system, we 
can describe the system in terms of density n(x,t) of molecules at point x G [0, L] and 
time t, so that n(x, t) Sx gives the number of molecules in the small interval [x, x + 5x] 
at time t. The evolution of n(x, t) is governed by the diffusion equation 

Yt= D % forxe[0,L], t>0, (2) 

where D is the diffusion constant. The general first-order reactive boundary condition 
at x = is the so-called Robin boundary condition 

dn 

D—{0,t) = Kn{0,t) (3) 

where the constant K describes the reactivity of the boundary (see Appendix for the 
relation between K and the chemical properties of the boundary) and may in general 
depend on time. The boundary condition at right boundary x = L can be treated 
similarly. 

In the following four subsections we introduce four stochastic models of diffusion. 
The first model, introduced in Section [2TTI is a position jump process on a lattice. This 
model is discrete in both time and space, and is used in a stochastic simulation algorithm 
which is based on the react ion- diffusion master equation [101 [TTJ . The second model, 
introduced in Section 12.21 is again discrete in time and discontinuous in space but the 
positions of molecules are not confined to a regular lattice. It is essentially the Euler 
scheme for the Smoluchowski stochastic differential equation, which is the core of the 
stochastic approach of Andrews and Bray [lj. The third scheme, introduced in Section 
12. 3[ is a discrete velocity jump process which is a discrete in time, continuous in space 
random walk with discretized velocities, where the velocities evolve on a finite lattice. 
The last stochastic model of diffusion, introduced in Section I2.4[ is the Euler scheme for 
the solution of the stochastic Langevin equation. It is a velocity jump process (i.e. a 
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random walk discrete in time, continuous in space and discontinuous in velocities) where 
the Brownian particle can move with any real value of the velocity. In all four cases, we 
study the connections between the boundary conditions of the stochastic simulation and 
Robin boundary condition (j3J) of the macroscopic diffusion equation (jSJ). We provide 
the relation between K in ([3]) and the parameters of each model. The mathematical 
derivation of these relations is included later, in Section HI 

2.1. Position Jump Process I 

Let us discretize the domain of interest [0, L] into M lattice points a distance h = L/M 
apart, namely we consider the lattice 

(h 3h 5h 7h 9h (2M - l)h 

\2' ~2~' Y'T'T"" 2 

Let us choose time step At such that 2D At <C h 2 . We simulate a system of N molecules 
whose positions are assumed to be at one of the discrete positions (TjJ. Let Xi(t) be 
the position of the i-th molecule, i = 1, 2, . . . , N, at time t. The position Xi(t + At) is 
computed from the position Xi(t) as follows: 

(Xi(t) with probability 1 - 2DAt/h 2 , 

Xiit) - h with probability DAt/h 2 , (5) 
Xiit) + h with probability DAt/h 2 . 

At x = 0, we implement the following boundary condition: whenever a molecule hits 
the boundary, it is adsorbed with probability P\h, and reflected otherwise. Here, Pi is a 
given nonnegative constant. The implementation of this boundary condition at x = 
is performed as follows. If the i-th molecule is at position Xi(t) = h/2 and attempts to 
jump to the left, then 

x^t + At) = h/2 with probability 1 - P x h. (6) 

Otherwise, we remove the molecule from the system. We show in Section 14.11 that the 
random walk (jSJ) with boundary condition leads to the diffusion equation (j2J) with 
Robin boundary condition Q, where 

K 

K = PiD, which is equivalent to Pi = — . (7) 

2.2. Position Jump Process II 

Let us choose a time step At. Let X{(t), i = 1,2, ... ,N, be the position of the i-th 
molecule at time t. The position Xi(t + At) is computed from the position Xi(t) as 
follows: 

x^t + At) = Xi{t) + V2D At 0, i = l,...,N, (8) 

where Q is normally distributed random variable with zero mean and unit variance. This 
random walk is essentially the Euler scheme for the Smoluchowski stochastic differential 
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equation (|2~T1) as discussed later, in Section 14.21 We implement the following partially 
adsorbing boundary condition at x = 0: whenever a molecule hits the boundary, it 
is adsorbed with probability P 2 \fAt, and reflected otherwise. Obviously, if Xi(t + At) 
computed by (|HJ) is negative, a molecule has hit the boundary. However, Andrews and 
Bray pQ argue that there is a chance that a molecule hit the boundary during the finite 
time step even if x-i{t + At) computed by (jSJ) is positive that is, during the time interval 
[t, t + At] the molecule might have crossed to Xi negative and then crossed back to 
Xi positive again. They found that the probability that the molecule hit the boundary 
x = at least once during the time step At is exp[— Xi(t)xi(t + At)/(DAt)] for Xi(t) > 0, 
Xi(t+At) > 0. Consequently, the partially reflective boundary condition is implemented 
as follows: 



(a) If Xi(t + At) computed by (jHJ) is negative then Xi(t + At) = —Xi(t) — \JlD At Q 
with probability 1 — P 2 \f~At, otherwise we remove the molecule from the system. 

(b) If Xi(t + At) computed by (jSJ) is positive then we remove the molecule from 
the system with probability exp[— Xi(t)xi(t + At) / (DAt)}P 2 \fAt. 

The partially adsorbing boundary condition (a) - (b) leads to the Robin boundary 
condition ([3]) with 

2P 2 VD , , , KJnr 

K = — — , which is equivalent to P 2 = — 7=. (9) 



Let us note that some authors use the case (a) only as the implementation of the 
partially reflective boundary condition [20] . i.e. they do not take the Andrews and 
Bray correction (b) into account. Considering the random walk (jSJ) with the boundary 
condition (a) only, the parameter K of Robin boundary condition (|2D can be computed 

as 



K = 2 , which is equivalent to P 2 = — (10) 



Comparing ([9]) and (flOl) . we see that we lose a factor of 2 if we do not consider the 
Andrews and Bray correction (b). The mathematical justification of formulas (Q and 
(TTU|) is presented in Section 14.21 

2.3. Velocity Jump Process I 

We consider that each molecule moves along the x-axis at a constant (large) speed s, but 
that at random instants of time it reverses its direction according to a Poisson process 
with the turning frequency 

s 2 

x =h- <"> 

Therefore, the z-th molecule is described by two variables: its position Xi(t) and its 
velocity Vi(t) = ±s. We use a small time step At such that A At 1. During each time 
step a molecule moves with speed s in the chosen direction. At the end of each time 
step, uniformly distributed random number G [0, 1] is generated. If < X5t, then 
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the z-th molecule changes its direction, so that it will move during the next time step 
in the opposite direction. 

We implement the following partially adsorbing boundary condition at x = 0: 
whenever a molecule hits the boundary, it is adsorbed with probability P^/s, and reflected 
otherwise. Here, P3 is a given nonnegative number. The partially adsorbing boundary 
condition at x = is implemented as follows. If the position of the z-th molecule satisfies 
Xi(t + At) < at the end of the time step, then 

*f + *\ = -^7f A * } with probability!-^, (12) 
Vi(t + At) = - Vi (t) J P s V ; 

otherwise, the z-th molecule is removed from the system. It can be shown that this 
velocity jump process leads to diffusion equation (jSJ) provided that s is sufficiently 
large (see Section fl~3l for details). Boundary condition (112p can be related with Robin 
boundary condition (j3J), with 
P 3 

K = — , which is equivalent to P3 = 2K. (13) 
2.4- Velocity Jump Process II 

Let us choose time step At. The z-th molecule is described by two variables: its position 
Xi(t) and its velocity Vi(t). We compute position Xi(t + At) and velocity Vi(t + At) from 
position Xi(t) and velocity Vi(t) by formulas 

Xi(t + At) = Xi(t) +v t (t)At, (14) 
Vi(t + At) = v t (t) - f3vi(t)At + (3V2DAt(i, (15) 

where (3 is the (large) friction coefficient and Q is a normally distributed random variable 
with zero mean and unit variance. This random walk is essentially the Euler scheme 
for the stochastic Langevin equation [3] . The partially reflective boundary condition at 
x = can be stated as follows: whenever a molecule hits the boundary, it is adsorbed with 
probability P4/ y/]3, and reflected otherwise. Here, P4 is a given nonnegative number. The 
implementation of this boundary condition is straightforward. If Xi(t + At) computed 
by (jl4jl is negative then 

Xi (t + At) = - Xi (t) -Vi(t)At \ P 4 f . 

Vl (t + At) = -v t (t) + ^t)At-pV2DAtQ ] -thprobabrhtyl-^; (16) 

otherwise, we remove the z-th molecule from the system. It can be shown that this 
velocity jump process leads to diffusion equation (j2J) provided that (3 is sufficiently large 
(see Section for details). The parameter K of Robin boundary condition (J3]) is 



K = A , which is equivalent to P4 = — (17) 

V 2n V D 
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3. Comparison of stochastic models of diffusion 

In this section, we present the results of two illustrative numerical simulations. First, 
we choose the macroscopic value of K in ([3]), and we show the results of stochastic 
simulations with the correct choice of the probabilities of the adsorption on the reactive 
boundary for each model. We demonstrate numerically the validity of relations between 
these probabilities and K, which were stated in Section [2] (the mathematical justification 
of these formulas is provided later, in Section HJ). In the second numerical example, we 
choose the apparently same microscopic boundary condition for each model. The goal 
is to demonstrate that the realistic boundary condition has to be chosen for each model 
with care, by applying formulas ((Tj), Q), (TT3"]) or (fTTj) . 

3.1. Stochastic simulation of models of diffusion 

Let us consider the computational domain [0,5], i.e. L = 5 in this section. We choose 
diffusion constant D — 1, and reactivity of the boundary x = as K = 2. We 
consider that the right boundary x = L is reflecting. Given an initial density profile 
n(x, 0), we can compute the density n(x,t) at time t > by solving diffusion equation 
(J2]) accompanied with Robin boundary condition (j3J) at x = and no-flux boundary 
condition at x = L. In this section, we show that comparable results can be obtained by 
all four stochastic models provided that we choose the boundary conditions accordingly. 

The key formulas were provided in the previous section. Given values of K and D, 
we can compute the adsorbing probabilities on the reactive boundary x = by formula 
(|7|) for Position Jump Process I, formula (02) for Position Jump Process II, formula (JT3l) 
for Velocity Jump Process I and formula (fTTl) for Velocity Jump Process II. We obtain 
appropriate values of constants Pi, P 2 , P3 and P4 which are used in the corresponding 
stochastic model. In our case K = 2 and D = 1, so that formulas (J7J), ([9]), ( |T3l) and ( fTTl) 
imply 

P x = 2, P 2 = ^ = 1.772, P 3 = 4, P 4 = = 5.013. (18) 

The results of stochastic simulations are presented in Figure [TJ The initial condition 
was chosen as follows: we start with 100, 000 molecules in domain [0, 5]. We put 75, 000 
molecules to position x = 1 and 25, 000 to position x = 2 at time t = 0. We plot the 
density profile at time t = 1 in Figure [TJ To do that, we divide the interval [0, 5] into 
50 bins of length 0.1 and we plot the number of molecules in each bin at time t = 1 
(histograms). We also plot the solution of diffusion equation ([2]) accompanied by Robin 
boundary condition ([3]) at x = and no-flux boundary condition at x — L. We see that 
all four models of diffusion give the same results provided that we choose the partially 
adsorbing probability accordingly. 

To compute the simulation results from Figure [TJ we also had to specify the 
additional parameters of the stochastic models. We used space step h — 0.1 and time 
step At = 10 -4 in Position Jump Process I. We used time step At = 10~ 4 in Position 
Jump Process II. We used speed s = 40 and time step At = 10~ 5 in Velocity Jump 
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Position Jump Process I 



Position Jump Process II 
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Figure 1. Stochastic simulations of four different diffusion models for K = 2 and 
D = 1 . Probabilities of adsorption at partially adsorbing boundary x = were computed 
for each model according to formulas Q, (| 1 3\ and Q_?7| ). 



Process I and we used friction coefficient (3 = 200 and time step At = 10 6 in Velocity 
Jump Process II. 

3.2. Consequences of the same probability of adsorption 

Let us now consider the following boundary condition: whenever a molecule hits the 
boundary, it is adsorbed with probability R, and reflected otherwise. We can easily 
modify the computer codes which were used to compute stochastic simulation results 
from Figure Q] to incorporate this boundary condition: whenever the i-th molecule hits 
the boundary, we generate a random number rj uniformly distributed in [0, 1]. If fj < R, 
we remove the i-th molecule from the system. It means that the adsorbing condition 
which was stated in terms of P\, P2, P3 and P4 is replaced by the same condition 
expressed in terms of R. Choosing R = 0.05 and the same parameters and initial 
condition as in Section 13.11 the stochastic simulation results at time t = 1 are shown 
in Figure [2] (histograms obtained by dividing domain [0, 5] into 50 bins and plotting 
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Position Jump Process II 
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Figure 2. Stochastic simulations of four different diffusion models for R = 0.05 
and D = 1 (histograms) . Solid curves show the solution of the diffusion equation (TJj) 
accompanied by no-flux boundary condition at x = 5 and Robin boundary condition (T3j) 
at x — where the values of K are computed according to formulas Q, (j 13$ and 

([13- 



the number of molecules in each bin). We clearly see that the results quantitatively 
differ. The reason is that the probability of adsorption scales with other parameters of 
simulations: namely with space step h for Position Jump Process I, with time step At 
for Position Jump Process II, with speed s for Velocity Jump Process I and with friction 
coefficient (3 with Velocity Jump Process II. The values of these scaling parameters were 
chosen as in Section [3TTI Namely, we used space step h = 0.1 in Position Jump Process 
I, time step At = 10 -4 in Position Jump Process II, speed s = 40 in Velocity Jump 
Process I and friction coefficient f3 = 200 in Velocity Jump Process II. Using these values, 
we can compute P\, P2, P3 and P\ which correspond to R = 0.05. Moreover, we can 
use formulas ([TJ) , (J9j) , (|T3|) and (TTTj) to compute the corresponding value of K in Robin 
boundary condition ([3]). We obtain K = 0.5 for Position Jump Process 1, K = 5.64 
for Position Jump Process II, K — 1 for Velocity Jump Process I and K = 0.28 for 
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Velocity Jump Process II. The solutions of diffusion equation (J2J) accompanied by no- 
flux boundary condition at x = 5 and Robin boundary condition (j3J) at a; = with 
the appropriate choice of K are plotted in Figure [2] for comparison as solid curves. We 
see that the Robin boundary condition (j3J) at x = with the appropriate choice of K 
gives the correct results when compared with stochastic simulations. Moreover, we also 
confirm that the same value of R leads to the different values of K. Hence the boundary 
condition cannot be formulated in terms of one probability R. It has to be appropriately 
scaled as shown in Section [2j 

To enable a direct comparison between models, we can slightly reformulate Position 
Jump Process I. The formulation from Section I2TT1 was chosen in a way which is used in 
the stochastic reaction-diffusion approaches which are based on the reaction-diffusion 
master equation [TO], (TT]. In particular, no relation between h and At was given and 
the probability of partial adsorption had to be scaled with h. One can also formulate 
the position jump process on lattice as follows: we choose time step At and space 
step h = \JlD At. At each time step, the molecule jumps to the left with probability 
1/2 and to the right otherwise. This random walk can be accompanied with partially 
adsorbing boundary condition: whenever a molecule hits the boundary, it is adsorbed 
with probability Pi\/A£, and reflected otherwise. This boundary condition leads to Robin 
boundary condition ([3]) with K given by 

V2 

Comparing this formula with (Q or ffl3|) . we see that the Robin boundary condition is 
different for Position Jump Process I and Position Jump Process II even if we scale the 
adsorption probability with the same factor yAt. 

The velocity jump processes can also be further compared. To do this, let us note 
that the speed s of a molecule can be estimated as ^JkT/m where k is the Boltzmann's 
constant, T is the absolute temperature and m is the mass of a molecule. In particular, 
we get the relation s = y/D(3 which can be used to scale the boundary condition of 
Velocity Jump Process I in terms of y/]3 instead of s. Consequently, we can relate P3 
to P4 by P3 = P^\f~D. However, using this relation in ({TBI , we obtain a different Robin 
boundary condition fl3]) than by using formula (TT7j) . 

To summarize this section, it is possible to reformulate Position Jump Process I to 
have the adsorption probability of both position jump processes scaled as Py/At. It is 
possible to relate s to (3 in Velocity Jump Process I to have the adsorption probability 
of both velocity jump processes scaled as Pj Then all four cases lead to Robin 
boundary condition ([3]) of the form 

K = aPVD (19) 

where a is model-dependent. Consequently, the boundary condition has to be chosen 
differently for each model to get the same value of K in Robin boundary condition ([3]). 
One has to use formulas © , © , (jTBl and (fT7|) as we showed in Section 13.11 
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4. Mathematical analysis of stochastic models of diffusion 

The goal of this section is to provide the justification of the results from Section El For 
each stochastic model, we show that the model leads to the diffusion equation (j2J) away 
from the boundary. Moreover, we derive the Robin boundary condition for each model. 



4-1. Position Jump Process I 

Let Pk{t) be the probability of finding a molecule at mesh point xj~ = (2k — l)h/2 where 
k = 1, 2, . . . , M. If k ^ 1, k ^ M, then p k satisfies 

. . / 2DAt\ , , DAt ( , , . .\ 

Pk (t + At) = (1 - -p- ) Pk(t) + (p k+ i(t) +p k -i(t)j 

which can be rewritten as 

Pk(t + At) -pk(t) _ Pk+i{t) +p k -i(t) - 2p k (t) 
At h 2 

Passing to the limit At — > 0, h — > 0, we obtain the diffusion equation d2J). The boundary 
condition at x = can be incorporated into the equation for p\(t) as 

. . / 2DAt\ , , DAt ( / \ / . x , ,\ 
Pl (t + At) = f 1 - —^-\ p x {t) + (p 2 (t) + (1 - P 1 h)p 1 (t)j (20) 

which can be rewritten as 

— Pl (t + At)- Pl (t) _ DVAl f p 2 (t)- P i(t) 
At At " ~h— { h PlMt) 



Passing to the limit At —>■ 0, h —>■ such that yAt/h is kept constant, we obtain ([7j) 



4-2. Position Jump Process II 

Position jump process © is a discretized version (Euler scheme) of the stochastic 
differential equation 

X(t + dt) = X(i) + V2DdW(dt) (21) 

where dW(dt) is the normal random variable with mean and variance dt (i.e. the 
propagator of the special Wiener process) and D is the macroscopic diffusion constant. 
The diffusion equation (j2J) is the Fokker-Planck equation corresponding to the stochastic 
differential equation ([211 and its derivation can be found in any standard textbook 
|17j . The derivation of Robin boundary condition ([9]) is more delicate and requires 
the application of asymptotic methods [2]. To do that, we consider the Position Jump 
Process II on the semiinfinite interval [0, 00) subject to the boundary condition (a)-(b) 
from Section I2~2l Let p& t = PAt(x,t) : [0, 00) x AtN — > [0, 00) be the probability 
density function of the discretized process ([8]) with the boundary condition (a)-(b), so 
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that pAt(x,iAt)dx is the probability of finding a molecule in the interval [x,x + dx] at 
time t = iAt. We have 



poo 

p At (x,t + At) = / p(x, t + At | y, t) PAt(y, t) dy, 
Jo 



(22) 



where p(x, t + At \ y, t) is the conditional probability distribution function of finding 
a molecule at point x at time t + At given that it is at point y at time t. There 
are two possible options to reach point x at time t + At: either we use (JSj) only, 
i.e. x = y + \/2D At Q; or we use the boundary condition x = —y — \J2D At Q 
with probability 1 — P 2 V / At. If the former is true, we have to take into account that 
some molecules are lost because of the Andrews and Bray boundary correction (b). 
Consequently, we have 



p(x, t + At | y, t) 



1 



1 — exp 



VAnDAi 
+ (1 -P 2 v / At)exp 
1 



xy 



DAt 

jx + yf 
AD At 



P 2 v / Ai| exp 



{x - yf 
AD At 



VAnD At 
Thus ff22|) reads as follows 

PAt(x,t + At) = 

PAt(y,t) 
VAnD At 





(x-y) 2 ' 


^exp 


AD At _ 



[1 - 2P 2 VAt) exp 



{x + yf 

AD At 





{x-y) 2 ' 


^exp 


AD At _ 



+ (1 -2-PaVAt) exp 



[x + yf 
AD At 



dy. 



(23) 



Away from the boundary, a steepest descent approximation to the integral as At 
leads to the diffusion equation (T5]). However, as observed in [20J, in the vicinity of 
the boundary, such a calculation needs to be modified: there is a boundary layer of 
width v^At, and it is the solution in the boundary layer which determines the boundary 
condition of the diffusion equation. In the boundary layer, we change variables from x 
to 7] by setting x = a/AT r) and define the inner solution 

Pinner{v,t) = PAt(VAt 1],t) . 

Expanding Pi nner in the powers of V~At, we obtain 

Pinner 

(r), t) ~ p ifi (r], t) + VAtp iA (ri, t) + Atp ij2 (r], t) + . . . . 

Using this expansion in the integral equation (123]) and comparing the terms of the same 
order, we obtain at 0(1) that p i)0 is independent of r\. At 0(y/At) we find that 



+ 



PiAv) = 

° pm(Q 

VAnD 



-IP, 



exp 



Pi,o 
VAttD 

AD 



exp 



+ exp 



AD 

(Z + v) 2 

AD 



(24) 
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Now by matching the inner boundary layer expansion with the outer expansion 
PAt{x, t) ~ n(x, t) + . . ., we find Pi ; o(t) = n(0, t) and 



lim ?m m = p M . 

v ^oo or] ox 



(25) 



Thus to determine the boundary condition on n at x = we need to determine the 
behaviour of dpi t i/drj as i] — > oo. Differentiating (|24l) with respect to rj, we obtain 



dp 



v.i 



+ 



drj 



(v) 



PiV. 



'txD 



exp 



2 1 



V 
AD 



VAt^D V ZD 



V 



exp 



(g - vf 

AD 



i + V 
2D 



exp 



(Z + vY 

AD 



d£. 



Using integration by parts, we obtain the integral equation 



dpi, 



(V) 



P2 Pi,a 



+ 



1 



exp 



V 

AD 



(26) 



dp, 



drj 



VAnD Jo 
Let us define the function g(rj) by 

P% Pi,0 



(0 exp 









AD 


— exp 


AD 



exp 



r/- 



1 d P^i \ 



Then ( 1261) can be rewritten as 
0(77) = 0(r?) 
where 



VAnD 

<f>(v) 



9(0 (exp 

Pi Pi,o 

exp 



drj 

(g ~ ^) 2 
4L> 



exp 



4L> 



d£ 



2 i 



V8^D 

and the error function is defined by 



V 
8D 



erf 



erf 



'8D 



(27) 



(28) 



(29) 



2 p 

erf(f ) = — / exp [-a 2 ] da. 
V 71 Jo 



The function g(rj) is defined for 77 > 0. Since <p(rf) is odd function, we can define g(rj) for 
the negative values as an odd function too by setting g(j]) = —g(—7)) for r\ < 0. Then 
equation ( 1281) can be simplified to 



g (n) = 4>(r]) + 



VAnD 



9(0 exp 



AD 



de- 



(30) 



The natural way to solve such an equation is to apply a Fourier transform, but we have 
to be slightly careful since the Fourier transform of g does not exist in the classical sense 
(g tends to a constant at infinity, so is not integrable). Defining 

9+(v) = g(v)xio,oo){v) g-(v) = g(v)x(-oc,o](v), 
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where X[a,b] is the characteristic function of the interval [a, b] (that is, X[a,b](v) = 1 if 
a < 7] < b and zero otherwise), and applying the Fourier transform 

/oo 
^(77) exp[ikrj] dr/ 
-00 

to equation ( 1301) . we obtain 

gz{k) + r (fc) = m + {gz(k) + r m exp [-m 2 ] . 

Thus 

gW + ^)= l-exph^T <31) 
This seems like one equation for the two unknowns g+(k) and g~2(k), but in fact we know 
from their definitions that g+(k) is analytic for the imaginary part of k positive, while 
g~^{k) is analytic for the imaginary part of k negative, and this tells us how to divide 
all the poles of the right-hand side between g+(k) and g~^(k), except for the pole at the 
origin, which may appear in either g~+{k) or g~2(k). To divide this pole contribution up 
we note that since g is odd, g~+{k) = —g~2(—k), which implies that the coefficients of the 
odd powers of k near zero are equal in ~g+{k) and c/2(k), and that the coefficients of the 
even powers are zero. 
Using ( 1271) we have 

dpi 1 1 f°° ^ 

lim ' (77) = lim girf) = lim — / g~+{k) exp[— ikrf\ dk, 

where the inversion contour lies in the upper half-plane. Deforming the contour to — ioo 
we pick up residue contributions from each of the poles of (|3T!) in the lower half plane. 
The only finite contribution as rj — > 00 arises from the pole at the origin. Since 



we have 



6(h) , ±iVDP2Pi,ok 

'IX 



9+{k) ~ —7=^ as A; ^ 0, 
VvriJ k 

so that 

dp hl 2P 2 VD 
D hm -T—iv) = Pi,o- 

Using the matching condition (|25|) we therefore derive the Robin boundary condition 
(J9]) for the stochastic boundary condition (a)-(b). If we consider the boundary condition 
(a) only, equation ( |22l) leads to the modified formula (1231) where 2P 2 is replaced by P2. 
Thus the boundary layer method presented above gives in that case the Robin boundary 
condition ( flOi) . which differs from ([9]) by the factor of two. 
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4-3. Velocity Jump Process I 

Using standard methods jUJ [7], one can show that the density of molecules n(x, t) 
satisfies the damped wave (telegrapher's) equation 

1 d 2 n dn s 2 d 2 n , , 

+ 1^7 = TTT^- (32) 



2A dt 2 dt 2A dx 2 
The long time behaviour of fl32|) is described by the corresponding parabolic limit [22J. 
Using pij) . we obtain that n(x,t) satisfies (j2J) for times t ^> A -1 . 

Next, we derive the Robin boundary condition corresponding to (fl2l) . Let p + (x, t) 
be the density of molecules that are at (x, t) and are moving to the right, and let p~(x, t) 
be the density of molecules that are at (x, t) and are moving to the left. Then the density 
of molecules at (x,t) is given by the sum n(x,t) = p + (x,t) + p~(x,t), and the flux is 
j(x,t) = s(p + (x,t) — p~(x, t)). The stochastic boundary condition ( fl2i) implies 

P +(o,t)= (i- p ~(°^- 

This boundary condition can be written in terms of n and j as 

;( w+ iM).(,5);^.ia 

which implies 



Since 



P 3 n(0,t)=(^-2) j(0,t). 



j(0,t) n -D—(0,t), (33) 



dx 

we derive ([TBI) in the limit s — ► oo. 



Velocity Jump Process II 

The random walk (JT^ - liTBl) is a discretized version of Langevin's equation [3]. To derive 
the Robin boundary condition, we consider the behaviour of molecules in the semiinfinite 
interval [0, oo). The i-th molecule is described by two variables: its position Xi(t) and 
its velocity Vi(t). We compute the position Xi(t + At) and velocity Vi(t + At) from the 
position Xi(t) and velocity Vi(t) by (|T4|) - (fl5|) together with boundary condition (fl6|) at 
x = 0. Let f(x,v,t) be the density of molecules which are at position x with velocity 
v at time t, so that f(x,v,t)5x5v is number of molecules in interval [x,x + Sx] with 
velocities between v and v + 8v at time £. Assuming that the change in velocity of the 
i-th molecule during the time step is Av (i.e. Av = Vi(t + At) — Vi(t)), there are two 
possible options for the molecule to reach a point x > with velocity v at time t + At: 
either the molecule was at position Xi(t) = x — (v — Av)At with velocity Vi(t) = v — Av 
at time £; or at position Xi(t) = —x + (v + Av)At with velocity Vi(t) = —v — Av and 
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was reflected according to ([TBI) . Both cases make sense only if Xi(t) > 0. Consequently, 
f(x, v , t + At) can be computed from /(■, •, t) by the integral equation 



poo 

f(x,v,t + At) = f(x — (v — Av)At, v - Av, t) i>(v - Av; Av) dAv + (34) 

Jv—x/At 

+ (l - -^L J f f(-x + (v + Av) At, -v - Av, t)ip(-v - Av;Av) dAv 

where ip(w; Av) is a distribution function of the conditional probability that the change 
in velocity during the time step is Av provided that Vi(t) = w. Using (Tl5l) . we obtain 

(Av + (3wAtf 



ib(w: Av) = ; exp 

Py/AnDAt 



(35) 



A(3 2 DAt 

Passing to the limit At — > 0, we obtain that f(x,v,t) satisfies the Fokker-Planck 
equation [3] 



dt dx dv \ dv r 

together with the boundary condition 

f(0,v,t)= ^1-^/(0, -v,t). (37) 

The density of molecules at the point x and time t is defined by 

n(x,t) = I f(x,v,t)dv. (38) 
Jr 

To derive the diffusion equation for n and the corresponding Robin boundary condition 
we consider the limit in which (3 — > oo and rescale the velocity variable by setting 



v = riVPi f(x,T],t) = f(x,v,t), 

to give 

ldf 1 dj 9/ ? df\ 



(5 dt y 7 /? dx drj \ dr] 

We expand / in powers of 1/ \f]3 as 

- 11 

f{x,7],t) = fo{X,T],t) + -j=fx{x,T],t) +-f 2 (x,T),t) + .... (40) 

Substituting (j40l into fl39l) and equating coefficients of powers of (3 we obtain 

°U + D y>)=0, (41) 



drj \ dr\ 
drj \ drj J dx ' 

^^ /2+ VJ = ^ + ^- (43) 
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Solving equations (ll"Tl) - (l4"2j) . we obtain 



2D 



fo(x,r],t) = g(x,t)exp 
fi(x,r],t) = -^(x,t)r/exp 



7] 

2D 



(44) 
(45) 



where the function g(x,t) is independent of rj. Substituting ()44l) -(l45l) into ( 1431) gives 



d_ 

dr] 



Vf2 + D 



dr/ 



d 2 g 2 
W V 6XP 



V 

2D 



dg 

+ m exp 



rj- 
2D 



Integrating over rj gives the solvability condition 

dg = D d 2 g 
dt dx 2 

Using ( 1381) . we see that g(x,t) is proportional to density of individuals n(x,t) for large 
P. Consequently, n(x, t) satisfies the diffusion equation (j2J) for large f3. Multiplying (1571) 
by v and integrating over v, we obtain 

Pa 



VP 



vf(0, —v, t) dv 



where flux is defined by j(0,t) = Lvf(0,v,t)dv. Substituting 
derive the Robin boundary condition (fT 



into 



(46) 



we 



5. Boundary conditions for stochastic models of reaction-diffusion processes 

In this section, we show that reactions in the solution do not change the boundary 
conditions from Section [21 i.e. the boundary conditions of stochastic models of the 
react ion- diffusion processes are determined by the corresponding diffusion model. First, 
we illustrate this fact numerically in Sections 15.11 and 15.21 In Section 15.11 we use the 
stochastic approach based on the reaction-diffusion master equation [10l[TT], so that the 
corresponding diffusion model is the Position Jump Process I. In Section 15. 2\ we use 
the stochastic approach of Andrews and Bray [I], so that the corresponding diffusion 
model is the Position Jump Process II. Then, in Section 15. 3[ we provide mathematical 
justification of the fact that the presence of reactions in the solution does not influence 
the choice of the boundary condition. 



5.1. Nonlinear reaction kinetics 

We consider two chemicals A and B which diffuse in the domain of interest [0, 1] with 
diffusion constants Da and Db, respectively, and which react according to Schnakenberg 
reaction kinetics [18] - The chemical A is produced with a constant rate (from a suitable 
reactant which is supposed to be in excess) and degraded. The chemical B is also 
produced with a constant rate. Moreover, A and B react according to the cubic reaction 

2A + B 3A (47) 
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position position 

Figure 3. Stochastic simulations of the reaction- diffusion model of the Schnakenberg 
kinetics with the partially adsorbing boundary at x = (histograms). Panel on the left 
shows chemical A and panel on the right chemical B. Solution of ( |^i§| )-( pff| ) with the 
Robin boundary condition (|50p at x — and no-flux boundary condition at x = 1 is 
plotted as the solid line. 



where k c is reaction constant. We use a partially adsorbing boundary condition at x = 
and a reflective boundary condition at x — 1. 

The stochastic simulation algorithm is based on the reaction-diffusion master 
equation [10J E] . We divide the domain of interest into M compartments of the length 
h = 1/M which are assumed to be well-mixed. In particular, one can use the classical 
Gillespie's algorithm [9] to simulate stochastically the reactions in each compartment. 
The system is then described by two M-dimensional vectors [Ai, A 2 , . ■ ■ , Am] and 
[Bi,B 2 , . . . ,Bm] where Ai (resp. Bi) denotes the number of molecules of chemical A 
(resp. B) in the i-th compartment. The diffusion of chemicals is added to the system as 
another set of reactions-jumps between the neighbouring compartments with the rate 
Da/Ii 2 (resp. D B /h 2 ) [15j. In particular, the model of diffusion is equivalent to the 
Position Jump Process I. 

We choose M = 50 in what follows, i.e. h = 0.02. Initially, we have u = 1000 
molecules of A and B in each compartment, i.e. A^ = B^ = u, i = 1,2, ... , M, at 
time t = 0. The rate of production of A is 2u molecules per compartment per unit of 
time. The rate of production of B is 8a; molecules per compartment per unit of time. 
The degradation rate of A in the i-th compartment is proportional to 6Ai and k c is 
chosen to be 3uj~ 2 . Diffusion constants are Da = 1 and Db = 0.1. We implement 
the following boundary condition at x = 0: whenever a molecule of chemical A (resp. 
B) hits the boundary, it is adsorbed with probability P^h (resp. Pfh), and reflected 
otherwise. We choose = P^ = 10. We consider the reflective boundary condition 
for both chemicals at right boundary x = 1. Number of molecules in each compartment 
at time t = 1 are plotted in Figure [3] (histograms) . 
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Since oj is chosen sufficiently large, we can compare the results of stochastic 
simulation with A = ua and B = ub where a, b are the solution of the system of 
react ion- diffusion equations 

^ = D A ^ + 2-6a + 3an (48) 
at ox 2 

db d 2 b 

g^ + 8-3^ (49, 

with Robin boundary conditions at x = given by (j7j), namely 

^(0,f) = Pfa(0,t), f*(0,*) = Pf 6(0, t), (50) 
ox ox 

and with no-flux boundary conditions at right boundary x = 1. The curves A = cua and 

B = ub at time t — 1 are plotted in Figure [3] as solid lines for comparison. We see that 

the Robin boundary (j7]) which was derived for the corresponding diffusion model gives 

good results for the full react ion- diffusion simulation as well. 

We note that the system (j4"8]) - (j4"9"]) possesses a so-called Turing instability [13] if 

the values of the diffusion constants Da and Db are chosen to be sufficiently different. 

Our choice Da = 1 and Db = 0.1 falls in the regime in which the homogeneous solution 

a = 10/6 and b = 48/50 of ( l48l) -( l49l) is stable. On increasing the ratio Da/Db Turing 

patterns would develop, and we would observe solutions with multiple peaks (provided 

that the domain size is sufficiently large); see e.g. 



5.2. Spatially localized reactions 

In some morphogenesis applications [T9l [16] . one assumes that some prepatterning in the 
domain exists and one wants to validate the react ion- diffusion mechanism of the next 
stage of the patterning of the embryo. In this section, we present an example motivated 
by this approach. We consider one- dimensional domain [0, 3] and we use the molecular 
based approach of Andrews and Bray [Tj, i.e. the diffusion model is given by Position 
Jump Process II. We choose a small simulation time step At. At each time step, a 
molecule is released at random points in the subinterval [1, 2] with probability k p At -C 1. 
Moreover, we assume that any molecule is degraded with probability k^At <C 1 during 
the simulation time step. Here, k p and k^ are given constants. We implement the 
following boundary condition at x = 0: whenever a molecule hits the boundary, it is 
adsorbed with probability P2V / At, and reflected otherwise. We consider the reflective 
boundary condition at right boundary x = 3. 

We start with no molecules in computational domain [0, 3] at time t = 0. We choose 
diffusion constant D — 1, time step At = 10~ 7 , production rate k p = 10 5 , degradation 
rate kd = 1 and adsorption probability constant Pi = 5. To visualize the results, we 
divide the interval [0, 3] into 30 bins of length 0.1 and we plot the number of molecules in 
each bin at time t — 1 in Figure H] (histogram). The results of the stochastic simulation 
can be compared with the solution of reaction-diffusion equation 
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3500 




1 2 
position 



Figure 4. Stochastic simulations of the reaction- diffusion model with spatially localized 
reaction with the partially adsorbing boundary at x = (histogram). Solution of 
<\51\ with Robin boundary condition (0), with K = 5.642, and with no-flux boundary 
condition at x = 3, is plotted as the solid line. 



-qI = D Q^ + Qlk pX[x,-A -k d n, (51) 

where X[i,2] is characteristic function of the interval [1,2]. Equation (15 ip is solved in 
the interval [0,3] together with Robin boundary condition ([3]), with K given by (Q, 
and with no-flux boundary condition at x — 3. Using P 2 = 5, formula implies that 
K = 5.642. The density profile n(x, 1) at time t — 1 is plotted in Figure H] as a solid 
line for comparison. We see that the Robin boundary condition which was derived 
for the corresponding diffusion model, gives good results for the full reaction-diffusion 
simulation algorithm too. 



5. 3. Mathematical justification 

Let us consider a system of k chemicals diffusing and reacting in domain [0, L\. Let us 
suppose that the diffusion model is the Position Jump Process I, i.e. we use the stochastic 
approach based on the reaction-diffusion master equation [TO] |TT] as in Section 15. 11 Let 
p\{t) (resp. p 3 2 (t)) be the number of molecules of the j-th chemical, j = 1,2, ... ,k, at 
the boundary mesh point x\ = h/2 (resp. at x 2 = 3h/2). If there are no reactions going 
on, then p{(t) and p 2 {t) are related according to (|20|) . Introducing reactions at mesh 
point x = h/2, the boundary equation (1201 is modified as follows 

m + ^) = (i - ^) p{(t) + ^ (pat) + (i - p 1 hy l (t)) + Mf( P \(t), . . . ,p\{t)) 

where f(p\(t), . . . ,p\(t)) is the sum of the rates of all reactions which modifies the j-th 
chemical. Following the same procedure as in Section [4711 we find out that the additional 
term does not influence the boundary condition (it is O(At) and only 0(\ /r At) terms 
have nonzero contribution to the Robin boundary condition). In particular, we conclude 
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that the Robin boundary condition of the stochastic reaction-diffusion model is given 
by©- 

Next, let us consider the stochastic reaction-diffusion model of Andrews and Bray 
PQ which was used in Section l5\2l Let p 1 At = p J M {x,t) : [0, oo) x AtN — > [0, oo) be the 
density function of the j-th. chemical species, so that p& t (x, iAt)dx is the number of j-th 
molecules in interval [x, x + dx] at time t = iAt. If there are no reactions going on, then 
p J At satisfies the formula fl22|) . Introducing the reactions to the system, formula (}2"2"|) is 
modified as follows: 

POO 

p At (x,t + At) = [l + 0(At)}p(x,t + At\y,t)p At (y,t)dy (52) 
Jo 

where the additional O(At) term corresponds to the reactions in the solution. As before, 
this term is of lower order (compared to O(vAt) terms) and does not influence the Robin 
boundary condition. Consequently, the Robin boundary condition ([3]) is obtained with 
K given by ([9]). 

In this paper, we did not use the velocity jump processes to simulate stochas- 
tically the reaction-diffusion process. Velocity jump processes are generally more 
computationally intensive. However, they might be of use if one considers that only 
sufficiently fast molecules can actually react. Alternatively, one can use the approach 
based on binding/unbinding radii, as in the Andrews and Bray method [lj, to incorporate 
higher order reactions to the velocity jump models of molecular diffusion. In any 
case, the reactions are adding again O(At) terms and do not influence the boundary 
conditions, i.e. the boundary conditions of the stochastic react ion- diffusion models can 
be chosen as in the corresponding model of the diffusion only. 

6. Conclusions and outlook 

We have derived the correct boundary conditions for a number of stochastic models of 
react ion- diffusion processes. For each model, we related the (microscopic) probability of 
adsorption on the boundary with the (macroscopic) Robin (reactive) boundary condition 
(j3J). First, we studied several stochastic models of diffusion. We showed that each model 
is suitable for the description of the molecular diffusion far from the boundary. Moreover, 
we derived formulas (EI) , (J9j) , ( |T3i) and ( |T71) relating reactivity K of the boundary with the 
parameters of the stochastic simulation algorithms. Then, we showed that the boundary 
conditions of stochastic models of reaction-diffusion processes are the same as for the 
corresponding model of diffusion only. We studied the stochastic approaches based on 
the reaction-diffusion master equation [TUl [11] and on the Smoluchowski equation pp. 
The main conclusion of this work is that a modeller can use any of the stochastic model 
of the diffusion from Section [2, provided that the adsorption probability on the reactive 
boundary is chosen according to the corresponding formula, i.e. (0), (EJ), (fi~3l) or ( TT7T) . 
which is mo del- dependent. 

We also presented the mathematical derivation of key formulas ([7]), ([9]), ffT3l) and 
( fT71) . We devoted the most space to the derivation of formula ([9]) which is (to our 
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knowledge) a new mathematical result. Derivation of formulas fl7|) and (1131) is simple 
and we included the mathematical arguments for completeness. The last formula ffTTj) 
has already appeared in literature [Ti] . though our derivation is more systematic. 

It is interesting to note that in some applications, the reactivity of the boundary 
depends also on the geometrical constraints on the boundary. The binding sites 
on the surface (e.g. reactive groups or receptors) become full as the adsorption 
progresses. Moreover, attaching large molecules to a binding site can sterically shield 
the neighbourghing binding sites on the surface. For example, in (6j H] we studied 
the chemisorption of polymers where the attachment of a long polymer molecule to 
the surface prevents attachment of other reactive polymers next to it. This steric 
shielding was modelled using random sequential adsorption [8J. In these models, an 
adsorption of one molecule to the surface is attempted per unit of time. To relate the 
time scale of random sequential adsorption algorithms with physical time, one should 
couple the theory of reactive boundaries presented here with algorithms which take the 
additional geometrical constraints on the boundary into account. This is an area of 
ongoing research [5]. 
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-^7 = D—-Kn6(x), for x G [0, L], t > 0, (53) 



Appendix. Robin boundary condition and chemistry 

In this appendix, we show the relation between the Robin boundary condition ([3]) and 
the experimentally measurable chemical properties of the boundary. Let us consider a 
chemical diffusing in the domain [0, L] which is adsorbed by boundary x = with some 
rate K. This problem can be described by the reaction-diffusion equation 

dn d 2 n — 

dt dx 2 

together with no-flux boundary conditions, where D is the diffusion constant and 5(x) 
a Dirac delta function. The term Knd~(x) is a standard description of reaction kinetics, 
localized on the boundary. In the paper, we used an alternative description of the 
chemically adsorbing boundary, given by the diffusion equation (T5]) accompanied by the 
Robin boundary condition ([3]). It is interesting to note that the constant K in is 
actually equal to the experimentally measurable constant K. To see this, we discretize 
( |53l) with space step h and we denote n (t) = n{h/2,i) and ni(t) = n(3h/2,t). Using a 
no-flux boundary condition (i.e. n(—h/2,t) = n(h/2,t)), the discretization of (|53|) gives 

dn ni-n 1 , , 

~dt = h? KU0 h (54) 

which is equivalent to 

~dT = D h? • (55) 
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The same equation can be obtained by the discretization of the diffusion equation (}2"1) 
together with the Robin boundary condition ([3]). Hence we showed that K = K, i.e. 
the Robin boundary condition ([3]) is indeed the correct macroscopic description of the 
chemically reacting boundary. 
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